System identifying method

ABSTRACT

Quick real-time identification and estimation of a time-non-varying or time-varying system. A new H ∞  evaluation criterion is determined, a fast algorithm for a modified H ∞  filter based on the criterion is developed, and a quick time-varying system identifying method according to the fast H ∞  filtering algorithm is provided. By the fast H ∞  filtering algorithm, a time-varying system sharply varying can be traced with an amount of calculation O(N) per unit time step. The algorithm completely agrees with a fast Kalman filtering algorithm at the extreme of the upper limit value. If the estimate of impulse response is determined, a pseudo-echo is sequentially determined from the estimate and subtracted from the real echo to cancel the echo. Thus an echo chancellor is realized.

BACKGROUND OF THE INVENTION

[0001] The present invention relates to system identification methods, and more particularly, to a time-varying system identification method at a high speed in real time by using a fast algorithm for modified H_(∞) filters developed based on a new H_(∞) evaluation criterion.

[0002] In general, system identification is to estimate a mathematical model (a transfer function, an impulse response, or the like) of a system input-and-output relationship according to the input-and-output data. Typical applications thereof include echo cancellers in international communications, automatic equalizers in data communications, echo cancellers in acoustic systems, sound-field reproduction, and active noise control in vehicles. Details are written in “Digital Signal Processing Handbook”, the Institute of Electronics, Information and Communication Engineers, 1993, or the like.

[0003] (Basic Principle)

[0004]FIG. 14 shows a configuration for system identification. This system includes an unknown system 1 and an adaptive filter 2. The adaptive filter 2 has an FIR digital filter 3 and an adaptive algorithm 4.

[0005] One case which uses an output error method to identify the unknown system 1 will be described below. Here, u_(k) indicates an input to the unknown system 1, d_(k) indicates the output of the system, which is a signal to be obtained, and d{circumflex over ( )}_(k) indicates the output of the filter. (A mark “{circumflex over ( )}” means an estimated value and should be placed above characters, but it is placed at the upper right of the characters for input convenience. This notation may be used through the present specification.)

[0006] Since an impulse response is generally used as a parameter of an unknown system, the adaptive filter adjusts a coefficient of the FIR digital filter 3 by the adaptive algorithm such that an evaluation error e_(k)=d_(k)−d{circumflex over ( )}_(k) in the figure is minimized.

[0007]FIG. 15 shows the structure of an impulse-response adjustment mechanism.

[0008] Here, as an example of adaptive algorithm, the following LMS algorithm (least mean square algorithm) is widely used because of its computational simplicity.

[0009] [LMS Algorithm]

ĥ _(k+1) =ĥ _(k) +μu _(k)(y _(k) −u _(k) ^(T) ĥ _(k))  (1)

[0010] where,

ĥ _(k) =[ĥ ₀ [k], . . . , ĥ _(N−1) [k]] ^(T) , u _(k) =[u _(k) , . . . , u _(k−N+1)]^(T) , u>0  (2)

[0011] Generally, Kalman filters, which converges quickly, are suitable for identifying a time-varying system.

[0012] [Kalman Filter]

{circumflex over (x)} _(k|k) ={circumflex over (x)} _(k|k−1) +K _(k)(y_(k) −H _(k) {circumflex over (x)} _(k|k−1))

{circumflex over (x)} _(k+1|k) ={circumflex over (x)} _(k|k)  (3)

K _(k) ={circumflex over (P)} _(k|k−1) H _(k) ^(T)(1H _(k) {circumflex over (P)} _(k|k−1) H _(k) ^(T))⁻¹  (4)

{circumflex over (P)} _(k|k) ={circumflex over (P)} _(k|k−1) −K _(k) H _(k) {circumflex over (P)} _(k|k−1)

[0013] $\begin{matrix} {{\hat{P}}_{k + {1/k}} = {{\hat{P}}_{k/k} + {\frac{\sigma_{w}^{2}}{\sigma_{v}^{2}}I}}} & (5) \end{matrix}$

[0014] where,

{circumflex over (x)} _(k|k) =[ĥ ₀ [k], . . . , ĥ _(N−1) [k]] ^(T) , H _(k) =[u _(k−1) , . . . , u _(k−N)]

{circumflex over (x)} _(0|1)=0, {circumflex over (P)} _(0|−1)=ε₀ I, ε ₀>0  (6)

[0015] Here, the impulse response {h_(i)} of unknown system is obtained as a state estimate x{circumflex over ( )}_(k|k), and an input {u_(k)} to the system is used as an element of an observation matrix H_(k).

[0016] Fast Kalman filtering algorithm is also known, which reduces the amount of calculation per unit time step to the number of times calculations are performed proportional to N, that is, O(N) by applying the shift property (H_(k+1)(i+1)=H_(k)(i)) of the observation matrix H_(k) to a Kalman filter obtained when σ² _(w)=0. Details are written in “Digital Signal Processing Handbook”, the Institute of Electronics, Information and Communication Engineers, 1993, or the like.

[0017] (Applications to Echo Cancellers)

[0018] Four-wire circuits are used for long distance telephone lines such as for international calls for a reason of signal amplification and others. On the other hand, two-wire circuits are used for subscriber lines because they are relatively short.

[0019]FIG. 16 is an explanation view of a communication system and an echo. Impedance matching is performed by disposing hybrid transformers at connection points between two-wire circuits and a four-wire circuit, as shown in the figure. If the impedance matching is perfect, a signal (voice) from a speaker B reaches only a speaker A. However, it is difficult to make the matching perfect in general. A phenomenon occurs in which a part of a received signal leaks to the four-wire circuit, is amplified, and returns to the receiver (speaker A). This is an echo. As the transmission distance gets longer (delay time gets longer), the effect of the echo gets larger, and the quality of telephone speech significantly deteriorates (in pulse transmission, an echo influences significantly even in a short distance, and the quality of telephone speech deteriorates).

[0020]FIG. 17 shows a basic principle of an echo canceller.

[0021] As shown in the figure, the echo canceller is introduced to successively estimate the impulse response of an echo path by using a received signal and an echo directly observable, and to subtract a quasi echo obtained by using the estimate from an actual echo to cancel and remove the echo.

[0022] The impulse response of the echo path is estimated such that the mean square error of a residual echo e_(k) is minimized. Elements which disturb the estimation of the echo path are line noise and a signal (voice) from the speaker A. When two speakers start talking at the same time (double talking), the estimation of the impulse response is generally halted. Since the impulse response of the hybrid transformers has a length of approximately 50 [ms], if the sampling period is set to 125 [μs], the order of the impulse response of the echo path actually becomes approximately 400.

SUMMARY OF THE INVENTION

[0023] In conventional arts, the LMS algorithm (least-mean-square algorithm) has been widely used as adaptive algorithm due to its simplicity, but it is impossible to closely track a time-varying system which varies suddenly due to its very slow convergence.

[0024] A Kalman filter, which has an excellent tracking performance, requires the amount of calculation proportional to O(N²) or O(N³). Since the amount of calculation rapidly increases with the tap number N, it is difficult to process actual issues requiring a high tap number N in real time. As a countermeasure therefor, a fast Kalman filter with the computational complexity of O(N) per unit time step for a tap number N has been proposed, but it is impossible for the filter to identify a time-varying system because of its stationary characteristic (incapable of taking system noise into account).

[0025] In view of the above-described points, it is an object of the present invention to implement a fast real-time identification of time-varying and time-invariant systems by using a fast algorithm for modified H_(∞) filters developed based on a new H_(∞) evaluation criterion. It is another object of the present invention to include, as a particular case of the present algorithm, a fast Kalman filtering algorithm, and to determine theoretically the covariance of system noise which is dominant in the tracking performance of time-varying systems. It is still another object of the present invention to provide a fast time-varying system identification method which is substantially effective even when an input signal is discontinuously varied, such as in an echo canceller for a time-varying system which varies extremely as sudden line switching. It is a further object of the present invention to provide a system identification method which is applicable to echo cancellers in communication systems and acoustic systems, sound-field reproduction, and noise control.

[0026] In order to solve the problems described above, a new H_(∞) evaluation criterion is introduced, a fast algorithm for the modified H_(∞) filters is developed according to the criterion, and a fast time-varying system identification method based on the fast algorithm is proposed in the present invention. The fast algorithm according to the present invention is capable of tracking a time-varying system which varies rapidly, with the complexity of O(N) per unit time step. Further, it has a convenient property that it completely matches the fast Kalman filtering algorithm at a limit of γ_(f)=∞.

BRIEF DESCRIPTION OF THE DRAWINGS

[0027]FIG. 1 is a flowchart of a fast algorithm.

[0028]FIG. 2 is an explanation view (1) of the amount of calculation in each part of a modified H_(∞) filtering algorithm.

[0029]FIG. 3 is an explanation view (2) of the amount of calculation in each part of the modified H_(∞) filtering algorithm.

[0030]FIG. 4 is an explanation view of the amount of calculation when the order of matrix calculation is changed.

[0031]FIG. 5 is an explanation view (1) of the amount of calculation in a fast H_(∞) filtering algorithm.

[0032]FIG. 6 is an explanation view (2) of the amount of calculation in the fast H_(∞) filtering algorithm.

[0033]FIG. 7 is a view showing values of an impulse response {h_(i)} to be estimated.

[0034]FIG. 8 is a comparative explanation view of estimated results of impulse responses obtained by the modified H_(∞) filter and the fast H_(∞) filter.

[0035]FIG. 9 is a view of measurement results of computation time.

[0036]FIG. 10 is a view (1) of simulation results of each algorithm.

[0037]FIG. 11 is a view (2) of simulation results of each algorithm.

[0038]FIG. 12 is a view showing the relationship between γ, and ρ.

[0039]FIG. 13 is a view showing the relationship among tap numbers and computation times of impulse responses obtained by the fast H_(∞) filter, fast Kalman filter, and LMS algorithm.

[0040]FIG. 14 is a view showing a configuration for system identification.

[0041]FIG. 15 is a view showing the configuration of an impulse-response adjustment mechanism.

[0042]FIG. 16 is an explanation view of a communication system and an echo.

[0043]FIG. 17 is a view showing the principle of an echo canceller.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0044] Embodiments of the present invention will be described hereinafter. Details are shown, for example, in “Derivation of A Fast Algorithm of Modified H_(∞) Filters”, K. Nishiyama, IEEE international Conference on Industrial Electronics, Control and Instrumentation, RBC-II, pp.462-467, October, 2000.

[0045] 1. Description of Symbols

[0046] First, main symbols used in the embodiments of the present invention and whether they are known or unknown will be described.

[0047] x_(k): State vector or just state, unknown and to be estimated.

[0048] x₀: Initial state, unknown.

[0049] w_(k): System noise, unknown.

[0050] v_(k): Observation noise, unknown.

[0051] y_(k): Observation signal, known and input to a filter.

[0052] z_(k): Output signal, unknown.

[0053] H_(k): Observation matrix, known.

[0054] L_(k): Output matrix, known.

[0055] x{circumflex over ( )}_(k|k): State value of the state x_(k) at time k, estimated by using observation signals y₀-y_(k). Given by a filter equation.

[0056] x{circumflex over ( )}_(0|0): Initial estimate of a state, essentially unknown but set to 0 for convenience.

[0057] K_(s, k+1): Filter gain, obtained by matrix P{circumflex over ( )}_(k|k−1).

[0058] Σ_(wk): Corresponds to the covariance matrix of the system noise, known in theory but unknown in practice.

[0059] P{circumflex over ( )}_(k|k−1): Corresponds to the covariance matrix of the error of x{circumflex over ( )}_(k|k−1), given by a Riccati equation.

[0060] P{circumflex over ( )}_(1|0): Corresponds to the covariance matrix of an error in the initial state, essentially unknown but set to ε₀I for convenience.

[0061] σ² _(v): Variance of the observation noise, treated as known in theory but unknown in practice.

[0062] σ² _(w): Variance of the system noise, treated as known in theory but unknown in practice.

[0063] A mark “{circumflex over ( )}” placed above a symbol indicates an estimated value, a mark “U” indicates that the matrix is extended by one row, and a mark “{tilde over ()}” is added for convenience. These marks are placed at the upper right of characters for input convenience, but, as shown in expressions, they are identical with those placed above characters. “L”, “H”, “P”, and “K” indicate matrixes. Some of them are written in bold face as in expressions, but they are usually written in lightface for convenience.

[0064] 2. Modified H_(∞) Filter

[0065] Next, a state-space model as in the following equations (7) to (9) is discussed.

x _(k+1) =x _(k) +w _(k) , w _(k) ,x _(k) εR ^(N)  (7)

y _(k) =H _(k) x _(k) +v _(k) , y _(k) , v _(k) εR  (8)

z _(k) =H _(k) x _(k) , z _(k) εR, H _(k) εR ^(1×N)  (9)

[0066] where, L_(k)=H_(k) (H_(k)=[u_(k) u_(k−)2 . . . u_(k−N+1)]) assuming an echo canceller or the like. For such a state-space model, an H_(∞) evaluation criterion (γ_(f) is newly placed in the left-hand side) such as that shown by expression (10) is proposed. $\begin{matrix} {{\sup\limits_{x_{0},{\{ w_{i}\}},{\{ v_{i}\}}}\frac{\sum\limits_{i = 0}^{k}{{e_{f,i}}^{2}/\rho}}{{{x_{0} - {\overset{ˇ}{x}}_{0{- 1}}}}_{\underset{0}{\overset{- 1}{\Sigma}}}^{2} + {\sum\limits_{i = 0}^{k}{w_{i}}_{\underset{w_{k}}{\Sigma}}^{2}} + {\sum\limits_{i = 0}^{k}{{v_{i}}^{2}/\rho}}}} < \gamma_{f}^{2}} & (10) \end{matrix}$

[0067] When it is assumed that ρ or Σ_(wk) does not depend on γ_(f), a modified H_(∞) filter of level γ_(f) satisfying the evaluation criterion can be given by the following equations (11) to (14) by applying a standard H_(∞) estimation scheme known in the system identification field. This scheme is shown, for example, in “Linear Estimation in Krein Spaces—Part I: Theory,” B. Hassibi, A. H. Sayed, and T. Kailath, IEEE Trans. Automatic Control, 41, 1, pp.18-33, 1996., and “Linear Estimation in Krein Spaces—Part II: Applications,” B. Hassibi, A. H. Sayed, and T. Kailath, IEEE Trans. Automatic Control, 41, 1, pp.34-49, 1996. $\begin{matrix} {{\overset{\_}{z}}_{kk} = {H_{k}ϰ_{kk}}} & \quad & (11) \\ {{\hat{ϰ}}_{{k + 1}{k + 1}} = {{\hat{ϰ}}_{kk} + {K_{a,{k + 1}}\left( {y_{k + 1} - {H_{k + 1}{\hat{ϰ}}_{kk}}} \right)}}} & {{Filter}\quad {equation}} & (12) \\ {K_{a,{k + 1}} = {{\hat{P}}_{{k + 1}k}{H_{k + 1}^{T}\left( {{H_{k + 1}{\hat{P}}_{{k + 1}k}H_{k + 1}^{T}} + \rho} \right)}^{- 1}}} & {{Filter}{\quad \quad}{gain}} & (13) \\ {{{\hat{P}}_{{k + 1}k} = {{\hat{P}}_{k{k - 1}} - {{{\hat{P}}_{k{k - 1}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}{R_{e,k}^{- 1}\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}}{\hat{P}}_{k{k - 1}}} + {\sum_{w_{k}}{where}}}},} & {{Riccati}\quad {equation}} & (14) \\ \begin{matrix} {e_{j,i} = {z_{ii} - {H_{i}ϰ_{i}}}} \\ {R_{e,k} = {R_{k} + {\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}{P_{k{k - 1}}\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack}}}} \\ {{R_{k} = \begin{bmatrix} \rho & 0 \\ 0 & {- {\rho\gamma}_{j}^{2}} \end{bmatrix}},{\Sigma_{w_{k}} = {\gamma_{f}^{- 2}P_{{k + 1}k}}}} \\ {{{P_{k{k - 1}}^{- 1} + {H_{k}^{T}H_{k}}} < 0},\quad {P_{10} = {ɛ_{0}I}},\quad {ɛ_{0} > 0}} \\ {{{0 < \rho} = {{1 - \gamma_{f}^{- 2}} \leq 1}},\quad {\gamma_{f} > 1}} \end{matrix} & \quad & (15) \end{matrix}$

[0068] Since a weight ρ in the evaluation criterion depends on an upper limit γ_(f) determined in advance, the above algorithm is essentially different from that applied to normal H_(∞) filters. The present algorithm controls a maximum energy gain from disturbances (having the initial state x₀, the system noise {w_(i)}, and the observation noise {ν_(i)}) weighted by ρ to a filter error {e_(f,i)} so as to be smaller than γ_(f) ². Therefore, the present algorithm is a robust filtering algorithm against the disturbances. This property is reflected by the tracking characteristic of a time-varying system. When γ_(f)→∞ is satisfied, ρ=1 and Σ_(wk)=0. In this time, the modified H_(∞) filter becomes a normal H_(∞) filter.

[0069] The main load for calculating the modified H_(∞) filter rises during the update of P{circumflex over ( )}_(k+1|k)εR^(N×N), which requires the amount of calculation in proportion to N² or N³. That is, an arithmetic operation of O(N²) per unit time step is required. Here, a tap number N matches the dimension of the state vector x_(k). Therefore, as the dimension of x_(k) increases, the computation time required to perform the modified H_(∞) filter increases rapidly. In order to overcome the drawback, the introduction of a fast algorithm of the modified H_(∞) filter is needed.

[0070] 3. Fast H_(∞) Filtering Algorithm

[0071] The calculation of the Riccati equation (covariance equation of a state estimation error) shown in equation (14) is dominant in the computational complexity of the modified H_(∞) filter. Therefore, to process the modified H_(∞) filter at a high speed, if the filter gain of equation (13) is directly determined without using the Riccati equation, the computational burden can be significantly reduced.

[0072] Since it is difficult to derive a fast algorithm for directly obtaining a filter gain K_(s,k)εR^(N×1), however, evolving an algorithm for fast calculating a gain matrix defined as follows is examined.

K _(k) =P _(k) C _(k) ^(T) εR ^(N×2)  (16)

[0073] where, $\begin{matrix} \begin{matrix} {P_{k} = {\left\lbrack {_{k}^{T}\Omega_{k}_{k}} \right\rbrack^{- 1} = \left\lbrack {\sum\limits_{i = 1}^{k}{\rho^{k - i}C_{i}^{T}W_{i}C_{i}}} \right\rbrack^{- 1}}} \\ {{\Omega_{k} = \begin{bmatrix} {\rho\Omega}_{k - 1} & 0 \\ 0 & W_{k} \end{bmatrix}},{\Omega_{1} = W_{1}},{W_{i} = {{\rho \quad R_{i}^{- 1}} = {\begin{bmatrix} 1 & 0 \\ 0 & {- \gamma_{f}^{- 2}} \end{bmatrix} \in ^{2 \times 2}}}}} \\ {{_{k} = \begin{bmatrix} C_{1} \\ \vdots \\ C_{k} \end{bmatrix}},{C_{i} = {\begin{bmatrix} H_{i} \\ H_{i} \end{bmatrix} \in {^{2 \times N}.}}}} \end{matrix} & (17) \end{matrix}$

[0074] Here, the following lemmas are formed.

[0075] Lemma 1

[0076] A matrix P_(k) satisfies the Riccati equation of (14). Therefore, when a gain matrix K_(k) is obtained, the filter gain K_(s, k) is immediately obtained from the following lemma.

[0077] Lemma 2

[0078] The filter gain K_(s, k) of the modified H_(∞) filter is obtained by using the gain matrix K_(k) as shown below. In practice, the gain matrix K_(k) can be fast calculated by the recursive method in Lemma 3.

K _(s,k) =G _(k) ⁻¹ {overscore (K)} _(k) , G _(k) =ρ+γ _(f) ⁻² H _(k) {overscore (K)} _(k) εR  (18)

[0079] where,

{overscore (K)} _(k)(i)=ρK _(k)(i, 1), i=1,2, . . . , N  (19)

[0080] Lemma 3

[0081] The gain matrix K_(k) is updated as follows.

K _(k+1) =m _(k) −B _(k) F _(k) ⁻¹μ_(k) εR ^(N×2)  (20)

[0082] Here, m_(k)εR^(N×2) and μ_(k) εR^(1×2) are obtained by dividing a matrix of K^(U) _(k) =Q ^(U) _(k) ⁻¹C^(U) _(k) as shown below. $\begin{matrix} {\begin{bmatrix} m_{k} \\ \mu_{k} \end{bmatrix} = {\begin{bmatrix} 0 \\ K_{k} \end{bmatrix} + {\begin{bmatrix} S_{k}^{- 1} \\ {A_{k}S_{k}^{- 1}} \end{bmatrix}\left\lbrack {c_{k}^{T} + {A_{k}^{T}C_{k}^{T}}} \right\rbrack}}} & (21) \end{matrix}$

[0083] Auxiliary variables A_(k)εR^(N×1), S_(k)εR, and B_(k)F⁻¹ _(k)εR^(N×1) are obtained as well.

[0084] In conclusion, the fast H_(∞) filtering algorithm can be summarized as below.

[0085]FIG. 1 shows a flowchart of the fast algorithm, where L indicates a maximum data length.

[0086] [Step 0] Set initial conditions of a recursive expression as follows, where ε₀ is a substantially large positive constant. ${K_{0} = 0},{A_{- 1} = 0},{S_{- 1} = \frac{\rho}{ɛ_{0}}},{D_{- 1} = 0},{{\hat{x}}_{0|0} = 0}$

[0087] [Step 1] Compare time k with the maximum data length L. When the time k is larger than the maximum data length, terminate the processing. When the time k is equal to or smaller than the maximum data length, the processing proceeds to the next step (a conditional statement can be removed, if unnecessary).

[0088] [Step 2] Determine A_(k) and S_(k) recursively as follows.

{overscore (e)} _(k) =c _(k) +C _(k) A _(k−1) εR ^(2×1)

A _(k) =A _(k−1) −K _(k) W _(k){tilde over (e)}_(k) εR ^(N×1)

e _(k) =c _(k) +C _(k) A _(k) εR ^(2×1)

S _(k) =ρS _(k−1) +e _(k) ^(T) W _(k) {tilde over (e)} _(k) εR

[0089] [Step 3] Calculate K^(U) _(k) as follows. ${\overset{ˇ}{K}}_{k} = {\begin{bmatrix} {S_{k}^{- 1}e_{k}^{T}} \\ {K_{k} + {A_{k}S_{k}^{- 1}e_{k}^{T}}} \end{bmatrix} \in ^{{({N + 1})} \times 2}}$

[0090] [Step 4] Divide K^(U) _(k) as follows. ${{\overset{ˇ}{K}}_{k} = {{\begin{bmatrix} m_{k} \\ \mu_{k} \end{bmatrix}m_{k}} \in ^{N \times 2}}},{\mu_{k} \in ^{1 \times 2}}$

[0091] [Step 5] Determine D_(k), and obtain a gain matrix K_(s, k+1) from K_(k+1) as follows:

η_(k) =c _(k−N) C _(k+1) D _(k−1)

D _(k) =[D _(k−1) −m _(k) W _(k)η_(k)][1−μ_(k) W _(k) η _(k)]⁻¹

K _(k+1) =m _(k) −D _(k)μ_(k)

{overscore (K)} _(k+1)(i)=ρK _(k+1)(i, 1), i=1, . . . , N

K _(s,k+1) =G _(k+1) ⁻¹ {overscore (K)} _(k+1) , G _(k+1)=ρ+γ_(f) ⁻² H _(k+1) {overscore (K)} _(k+1)

[0092] where, η_(k)εR^(2×1), D_(k)εR^(N×1), K_(k+1)εR^(N×2), K_(s, k+1)εR^(N×1), 0<ρ=1−γ_(f) ^(−2≦1, γ) _(f)>1.

[0093] [Step 6] Update the filter equation of the H_(∞) filter as follows.

{circumflex over (x)} _(k+1|k−1) ={circumflex over (x)} _(k|k) +K _(a, k+1)(y _(k+1) −H _(k+1) {circumflex over (x)} _(k|k))

[0094] [Step7] Put the time k forward (k=k+1). The processing returns to Step 2, and the subsequent processes are repeated as long as the data exists.

[0095] Lemma 6 (Existence Condition Suitable for Fast Processing)

[0096] The existence of the fast H_(∞) filter can be checked with the computational complexity of O(N) by using the following existence condition.

[0097] [Existence Condition]

−e{circumflex over (Ξ)} _(i)+ργ_(f) ²>0, i=0, . . . , k  (22)

[0098] where, $\begin{matrix} {{\varrho = {1 - \gamma_{f}^{2}}},{{\hat{\Xi}}_{i}\frac{H_{i}{\overset{\sim}{K}}_{i}}{1 - {H_{i}{\overset{\sim}{K}}_{i}}}}} & (23) \end{matrix}$

[0099] 4. Computational Complexity of the Present Fast Algorithm

[0100] Next, how the computational complexity of the fast H_(∞) filtering algorithm decreases, as compared with the computational requirement of the modified H, filtering algorithm, will be discussed. Only the number of multiplications is used for evaluating the amount of calculation of an equation, and is calculated by the following method.

[0101] Number of multiplications when a J-by-K matrix is multiplied by a K-by-L matrix is J×K×L (times).

[0102] Here, when three or more matrixes or vectors are multiplied, they are calculated from the left unless a direction is specially shown in the figure.

[0103] (Computational Complexity of the Modified H_(∞) Filtering Algorithm)

[0104]FIGS. 2 and 3 are views showing of the amount of calculation of each part of the modified H_(∞) filtering algorithm, where N indicates a tap number. In FIG. 3(a), a calculation for obtaining R_(e, k) ⁻¹ from R_(e, k) is ignored. Similarly, in FIG. 2(a), a calculation for obtaining (H_(k+1)P{circumflex over ( )}_(k+1|k)H^(T) _(k+1)+1)⁻¹ from (H_(k+1)P{circumflex over ( )}_(k+1|k)H^(T) _(k+1)+1) is also ignored.

[0105] As shown in FIGS. 2(a), 3(a), and 3(b), the amount of calculation of each of K_(s|k+1), R_(e, k), and P{circumflex over ( )}_(k+1|k) is in proportion to the square of the tap number. Therefore, the amount of calculation of the entire modified H_(∞) filtering algorithm is O(N²) per unit time step.

[0106]FIG. 4 is a view showing the amount of calculation required when the order of matrix calculations is changed. More specifically, FIG. 4 shows the amount of calculation required when the order of matrix calculations in the following part is changed in the Riccati equation, compared with FIG. 3(b).

[0107] Since the amount of calculation of the above-described part is proportional to the cube of the tap number, the amount of calculation of P{circumflex over ( )}_(k+1|k) is also in proportion to the cube of the tap number. Accordingly, the amount of calculation of the entire H_(∞) filter increases from the square to the cube of the tap number.

[0108] Since either algorithm requires the amount of calculation proportional to the square or cube of the tap number, however, the computational burden for carrying out the filter increases significantly as the tap number increases. In fact, since a tap number used in the field of communication engineering, for example, is approximately 400, the practical use of the algorithm becomes very difficult.

[0109] (Computational Complexity of the Fast H_(∞) Filtering Algorithm)

[0110]FIGS. 5 and 6 are views showing the amount of calculation in the fast H_(∞) filtering algorithm. In the expression of K^(u) _(k) in FIG. 5(b), S_(k) ⁻¹ is obtained from S_(k), but the calculation thereof is ignored. Similarly, in the expression of D_(k) in FIG. 6(b), a calculation for obtaining [1−μ_(k)W_(k)η_(k)]⁻¹ from [1−μ_(k)W_(k)η_(k)] is also ignored.

[0111] The amount of calculation in the entire present fast algorithm is O(N) per unit time step according to FIGS. 5 and 6. Therefore, the amount of calculation in the fast H_(∞) filtering algorithm is in proportion to the tap number. In this case, the amount of calculation (the number of multiplications) for performing the fast H_(∞) filter once is 28N+16 per unit step, and is approximately double the amount (multiplication frequency) of calculation required for a fast Kalman filter, that is 12N+3.

[0112] As described above, although the computational complexity proportional to the square or cube of the tap number is required for the modified H_(∞) filtering algorithm, the computational complexity of the present fast algorithm is smaller and proportional to the tap number.

[0113] 5. Echo Canceller

[0114] The advantage of the present invention will be examined, with an echo canceller being taken as an example.

[0115] An observation value {y_(k)} of an echo {d_(k)} is expressed in the following expression by an (time-varying) impulse response ({h_(i)[k]} of an echo path, where it is considered that a received signal {u_(k)} is an input signal to the echo path: $\begin{matrix} {{y_{k} = {{d_{k} + v_{k}} = {{\sum\limits_{i = 0}^{N - 1}{{h_{i}\lbrack k\rbrack}u_{k - i}}} + v_{k}}}},\quad {k = 0},1,2,\ldots} & (24) \end{matrix}$

[0116] where, u_(k) and y_(k) indicate, respectively, the received signal and the echo at time t_(k) (=kT, T is a sampling period); v_(k) indicates circuit noise having zero mean at time t_(k); and h_(i)[k] (i=0, . . . , and N−1) is a time-varying impulse responses assuming a gradual change, and the tap number N thereof is known. Once estimated values {h{circumflex over ( )}_(i)[k]} of the impulse response are obtained, a quasi echo is generated as follows by using the estimated values. $\begin{matrix} {{{\hat{d}}_{k} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}},\quad {k = 0},1,2,\ldots} & (25) \end{matrix}$

[0117] Subtracting this from the echo (y_(k)−d{circumflex over ( )}_(k)≈0) , the echo is cancelled, where u_(k−1)=0 when k−i<0.

[0118] From the above description, the echo canceller problem is equivalent to successively estimating the impulse response {h_(i)[k]} of the echo path from the received signal {u_(k)} and echo {y_(k)}, both of which are directly observable.

[0119] In general, when the H_(∞) filter is applied to an echo canceller, equation (24) has to be expressed by a state-space model formed of a state equation and an observation equation. In this case, since the state vector to be obtained is the impulse response {h_(i)[k]}, allowing a state vector x_(k) to fluctuate with w_(k), the following state-space model can be constructed for the echo path.

{circumflex over (x)} _(k+1) ={circumflex over (x)} _(k+) w _(k) , {circumflex over (x)} _(w) _(k) εR ^(N)  (26)

y _(k) =H _(k) x _(k)+ν_(k) , y _(k),ν_(k) εR  (27)

z _(k) =H _(k) x _(k) , z _(k) ε R, H _(k) εR ^(1×N)  (28)

[0120] where,

x _(k) =[h ₀ [k], . . . , h _(N−1) [k]] ^(T) , w _(k) =[w _(k)(1), . . . , w_(k)(N)]^(T)

H _(k) =[u _(k) , . . . , u _(k−1) ], L _(k) =H _(k)

[0121] Modified H_(∞) filtering algorithm and fast H_(∞) filtering algorithm for such a state-space model are the same as those described above. While the impulse response is estimated, if the occurrence of a transmission signal is detected, the estimation is generally stopped in the meanwhile.

[0122] Thus, when an estimate {h{circumflex over ( )}_(i)[k]} of the impulse response is obtained, the quasi echo is successively obtained therefrom as follows. $\begin{matrix} {{\hat{d}}_{k} = {{H_{k}{\hat{x}}_{k|k}} = {\sum\limits_{i = 0}^{N - 1}{{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}}} & (29) \end{matrix}$

[0123] Therefore, subtracting this from an actual echo to cancel the echo, an echo canceller is implemented. Here, an estimate error, e_(k)=y_(k)−d{circumflex over ( )}_(k), is called a residual echo.

[0124] 6. Evaluation for Time-Invariant Impulse Response

[0125] (Evaluation of Estimation Accuracy)

[0126] A modified H_(∞) filter and a fast H_(∞) filter are evaluated by simulation in a case in which the impulse response of an echo path is time-invariant (h_(i)[k]=h_(i)) and the tap number N thereof is 24. $\begin{matrix} {y_{k} = {{\sum\limits_{i = 0}^{23}{h_{i}u_{k - i}}} + v_{k}}} & (30) \end{matrix}$

[0127]FIG. 7 is a view showing values of the impulse response {h_(i)} in this case. v_(k) is stationary Gaussian white noise having zero mean and variance σ_(v) ² of 1.0×10⁻⁶, and a sampling period T is set to 1.0 for convenience.

[0128] The received signal {u_(k)} is approximated by a quadratic AR model as shown below.

u _(k)=α₁ u _(k−1)+α₂ u _(k−2) +w _(k′)  (31)

[0129] where, α₁=0.7, α₂=0.1+, and w_(k)′ is stationary Gaussian white noise having zero mean and variance σ_(w′) ² of 0.04.

[0130] The modified H_(∞) filter and the fast H_(∞) filter will be compared.

[0131]FIG. 8 includes views showing estimated results of the impulse responses of the modified H_(∞) filter and the fast H_(∞) filter (initial value x{circumflex over ( )}_(0|0)=0, estimated value x{circumflex over ( )}_(100|100) at 100th step, ε_(o)=20). FIGS. 8(a) and (b) show estimated results of both filters when γ_(f)=10⁵, and FIGS. 8(c) and (d) show estimated results thereof when γ_(f)=2.0. From the figures, performance on the estimation accuracy of both filters is equal. In other words, speeding-up does not reduce the estimation accuracy. Note that, if γ_(f) is too small, the existence condition of the filters is not satisfied. When γ_(f)=1.0×10⁵, the results are substantially equal to that of a fast Kalman filter. Therefore, it is found that the fast H_(∞) filtering algorithm includes the fast Kalman filtering algorithm and its convergence rate can be accelerated by adjusting γ_(f).

[0132] (Evaluation of Computation Time)

[0133] Next, the computation time required for the modified H_(∞) filter and that for the fast H_(∞) filter are evaluated under conditions where the impulse response of the echo path is time-invariant and the tap number is increased to 24, 48, 96, 192, and 384. Since dispersion may occur in one measurement, the average of four measurements was used. The values shown in FIG. 7 are used as impulse responses {h_(i)} in simulation, and impulse responses {h_(i)} thereafter (24≦k<N) are set to 0. The filter calculation is performed up to step 100. The computation time was measured by a command “etime” of MATLAB on a workstation (sparc, 60 MHz, 32 MB).

[0134]FIG. 9 is a view showing measurement results of the computation time. In the Riccati equation, matrix calculation is performed for a modified H_(∞) filter (2) such that the amount of calculation is in proportion to the square of the tap number, and matrix calculation is performed for a modified H_(∞) filter (1) such that the amount of calculation is in proportion to the cube of the tap number (see FIG. 3(b) and FIG. 4). In modified H_(∞) filters, since the computational complexity is in proportion to the square or cube of the tap number depending on the order of matrix calculation as described above, they are not practical.

[0135] 7. Evaluation for Time-Varying Impulse Response

[0136] (Evaluation of Tracking Performance)

[0137] The tracking performance of each algorithm will be evaluated by using the echo canceller in a case in which the system (impulse response) is varied with time. It is assumed that the tap number of the impulse response is 48, and {h_(i)} is varied with time, as shown in FIG. 10(a), based on the values shown in FIG. 7. It is also assumed that v_(k) is stationary Gaussian white noise having zero mean and variance σ_(v) ² of 1.0×10⁻⁶, and the sampling period T is for convenience. The received signal {u_(k)} is approximated by a quadratic AR model as follows.

u _(k)=α₁ u _(k−1)+α₂ u _(k−2) +w _(k)′  (32)

[0138] Here, α₁=0.7, α₂=0.1, and w_(k)′ is stationary Gaussian white noise having zero mean and variance σ² _(w′) of 0.04.

[0139]FIGS. 10 and 11 are views showing the simulation result of each algorithm. These views show the tracking performance of time-varying systems which employ a fast H_(∞) filter (fast HF), a fast Kalman filter (fast KF), and LMS algorithm (LMS). FIG. 10(b) shows the estimates obtained with the fast H_(∞) filter when γ_(f)=2.0. FIG. 11(a) shows the estimates obtained with the fast Kalman filter. The initial value of the fast H_(∞) filter is set such that x{circumflex over ( )}_(0|0)=0 and so ε_(o)=20, and the initial value of the fast Kalman filter is set in the same way. FIG. 11(b) shows the estimates obtained by the LMS algorithm, wherein the initial value is set such that h{circumflex over ( )}_(o)=0, and the step size is set such that μ=0.5 so as to give a stable and rapid convergence. It is found that the tracking performance of the fast H_(∞) filter is extremely excellent, and the estimates become stable in about thirty steps after the impulse response is varied. On the other hand, the fast Kalman filter and the LMS algorithm cannot track the impulse response at all.

[0140] Generally, the tracking performance of H_(∞) filters having no system noise drops with time since the filter gain becomes smaller due to a decay in the diagonal component of P{circumflex over ( )}_(k|k−1) and the amount of update of the estimates decreases. In other words, as the number of steps increases, the estimates are updated little. Therefore, in order to improve the tracking performance of Kalman filters and H_(∞) filters, an appropriate value needs to be externally added to the diagonal component of the matrix P{circumflex over ( )}_(k|k−1). If it is directly added, however, a fast algorithm which uses the shift property of an observation matrix H_(k) cannot be implemented. It is one of significant features of the present invention to solve this problem theoretically by applying a weight ρ of 1−γ_(f) ⁻² to the H_(∞) evaluation criterion. The weight ρ appears in an update equation of S_(k) of the fast H_(∞) filtering algorithm, as follows.

[0141] (Update of Auxiliary Variable S_(k) of the Fast H_(∞) Filter)

[0142] An auxiliary variable S_(k) of the fast H_(∞) filter is indicated by the following expression.

S _(k) =ρS _(k−1) +e ^(T) _(k) W _(k) e{tilde over ()} _(k), 0<ρ=1−γ_(f) ⁻²≦1

[0143] In the fast H_(∞) filtering algorithm, S_(k) is used as S_(k) ⁻¹ in the equation of K^(u) _(k). In order to largely update the filter equation, S_(k) ⁻¹ must be larger. In other word, S_(k) needs to be kept small to make the large update. The existence of ρ prevents S_(k) from increasing rapidly, which is resultantly equivalent to adding system noise, and thereby the tracking performance is improved. Since the weight ρ is defined as 1−γ_(f) ⁻² the tracking performance can be varied by adjusting γ_(f) as confirmed in the simulation.

[0144]FIG. 12 is a view showing the relationship between γ_(f) and ρ. According to the figure, when γ_(f)=3.0, ρ=0.8889, which means that 89% of S_(k−)1 is transmitted to S_(k). Note that, if γ_(f) is set very small, however, the effect of S_(k−)1 is significantly reduced and the existence condition of the filter is not satisfied. When γ_(f) is large, γ=1. An increase in S_(k) is not suppressed at all, and therefore, the tracking performance drops. When γ_(f)=∞, in particular, the present fast algorithm completely matches the fast Kalman filtering algorithm.

[0145] (Evaluation of Computation Time)

[0146]FIG. 13 is a view showing the relationship among the tap number and the computation time for the fast H_(∞) filter, the fast Kalman filter, and the LMS algorithm, where the number of time steps executed for the filters is 300and γ_(f)=3.0. The computation time was measured for the fast H_(∞) filter, the fast Kalman filtering algorithm, and the LMS algorithm when the tap number was increased to 48, 96, 192, and 384 in the cases shown in FIGS. 10 and 11. Because dispersion may occur in one measurement result, the average of four measurement results, for example, was used.

[0147] It can be confirmed that, in any algorithm, the amount of calculation is in proportion to the tap number. It is also found that when the tap number is large, the computation time for the fast H_(∞) filtering algorithm is about a little less than twice the computation time for the fast Kalman filtering algorithm, and is approximately four times longer than that for the LMS algorithm, which is practical. Considering the tracking performance, it can be said that the fast H_(∞) filtering algorithm is sufficiently effective.

[0148] 8. Demonstration of Lemmas

[0149] Now, the above-described lemmas will be demonstrated.

[0150] (Demonstration of Lemma 1)

[0151] The inverse matrix of P_(k) will be indicated by equation (33). Further, a recursive equation for the matrix P_(k) can be obtained, as shown in equation (34), by using the matrix inversion lemma.

P _(k) ⁻¹ =ρO _(k−1) ^(T)Ω_(k−1) O _(k−1) +C _(k) ^(T) W _(k) C _(k)

=ρP _(k−1) ⁻¹ +C _(k) ^(T) W _(k) C _(k)  (33)

[0152] $\begin{matrix} {\begin{matrix} {P_{k} = \left\lbrack {{\rho \quad P_{k - 1}^{- 1}} + {\left\lbrack {H_{k}^{T}H_{k}^{T}} \right\rbrack {W_{k}\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}}}} \right\rbrack^{- 1}} \\ {= {{\rho^{- 1}P_{k - 1}} - {\rho^{- 1}{{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}} \cdot}}}} \\ {{{{\left( {W_{k}^{- 1} + {\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}\rho^{- 1}{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}}}} \right)^{- 1} \cdot \begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}}\rho^{- 1}P_{k - 1}},}} \end{matrix}\begin{matrix} {{\rho \quad P_{k}} = {P_{k - 1} - {{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}} \cdot}}} \\ {{{{\left( {R_{h} + {\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}}}} \right)^{- 1} \cdot \begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}}P_{k - 1}},}} \end{matrix}\begin{matrix} {P_{k} = {P_{k - 1} - {{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}} \cdot}}} \\ {{{{\left( {R_{h} + {\begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}}}} \right)^{- 1} \cdot \begin{bmatrix} H_{k} \\ H_{k} \end{bmatrix}}P_{k - 1}} + {\gamma_{f}^{- 2}{P_{k}.}}}} \end{matrix}} & (34) \end{matrix}$

[0153] It is understood, when P_(k) is replaced with P{circumflex over ( )}_(k+1|k), that the above equation satisfies the Riccati equation of (13).

[0154] (Demonstration of Lemma 2)

[0155] The gain matrix K_(k) can be expressed as follows. $\begin{matrix} \begin{matrix} {K_{k} = {{P_{k}C_{k}^{T}} = {\left\lbrack {{\rho \quad P_{k - 1}^{- 1}} + {C_{k}^{T}W_{k}C_{k}}} \right\rbrack^{- 1}C_{k}^{T}}}} \\ {{= {{\rho^{- 1}P_{k - 1}C_{k}^{T}} - {\rho^{- 1}P_{k - 1}{C_{k}^{T} \cdot}}}}\quad} \\ {{\left\lbrack {W_{k}^{- 1} + {C_{k\quad}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack^{- 1}C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \\ {= {{\rho^{- 1}P_{k - 1}C_{k}^{T}} -}} \\ {{\rho^{- 1}P_{k - 1}{{C_{k}^{T}\left\lbrack {W_{k}^{- 1} + {C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack}^{- 1} \cdot}}} \\ {\left\lbrack {\left( {W_{k}^{- 1} + {C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right) - W_{k}^{- 1}} \right\rbrack} \\ {= {\rho^{- 1}P_{k - 1}{C_{k}^{T}\left\lbrack {I + {W_{k}C_{k}\rho^{- 1}P_{k - 1}C_{k}^{T}}} \right\rbrack}^{- 1}}} \\ {= {\rho^{- 1}P_{k - 1}C_{k}^{T}{W_{k} \cdot \left\lbrack {W_{k} + {\rho^{- 1}W_{k}C_{k}P_{k - 1}C_{k}^{T}W_{k}}} \right\rbrack^{- 1}}}} \\ {= {\rho^{- 1}{{P_{k - 1}\left\lbrack {H_{k}^{T} - {\gamma_{f}^{- 2}H_{k}^{T}}} \right\rbrack}\left\lbrack {\begin{bmatrix} 1 & 0 \\ 0 & {- \gamma_{f}^{- 2}} \end{bmatrix} +} \right.}}} \\ \left. {{\rho^{- 1}\begin{bmatrix} H_{k} \\ {{- \gamma_{f}^{- 2}}H_{k}} \end{bmatrix}}{P_{k - 1}\begin{bmatrix} H_{k}^{T} & {{- \gamma_{f}^{- 2}}H_{k}^{T}} \end{bmatrix}}} \right\rbrack^{- 1} \\ {= {\rho^{- 1}{P_{k - 1}\begin{bmatrix} H_{k}^{T} & H_{k}^{T} \end{bmatrix}}\left( {1 + {H_{k}P_{k - 1}H_{k}^{T}}} \right)^{- 1}}} \end{matrix} & (35) \end{matrix}$

[0156] Further, the filter gain can be obtained from the first block column of the gain matrix K_(k), as shown in equation (18), by using G_(k)=(p+H_(k)P_(k−1)H^(T) _(k))/(1+H_(k)P_(k−1)H^(T) _(k)) and H_(k)K{tilde over ()}_(k)=H_(k)P_(k−1)H^(T) _(k)/(1+H_(k)P_(k−1)H^(T) _(k)).

[0157] (Demonstration of Lemma 3)

[0158] Assuming that the gain matrix K_(i) (i=0, . . . , and k) is given, the following matrix, K_(k+1), will be calculated.

Q _(k+1) K _(k+1) =C _(k+1) ^(T)  (36)

[0159] First, equations (37) and (38) are newly introduced to utilize the shift property of C_(k). Q^(U) _(k) is expressed recursively as shown in equation (39), and is divided as in the following equation (40). $\begin{matrix} {{\overset{˘}{C}}_{k}^{T} = {\begin{bmatrix} c_{k}^{T} \\ C_{k}^{T} \end{bmatrix} = {\begin{bmatrix} C_{k + 1}^{T} \\ c_{k - N}^{T} \end{bmatrix} \in ^{{({N + 1})} \times 2}}}} & (37) \\ {{\overset{˘}{Q}}_{k} = {{\sum\limits_{i = 1}^{k}{\rho^{k - i}{\overset{˘}{C}}_{i}^{T}W_{i}{\overset{˘}{C}}_{i}}} \in ^{{({N + 1})} \times {({N + 1})}}}} & (38) \end{matrix}$

 {haeck over (Q)} _(k) =ρ{haeck over (Q)} _(k−1) {haeck over (C)} _(k) ^(T) W _(k) {haeck over (C)} _(k)  (39 )

[0160] $\begin{matrix} {{\overset{˘}{Q}}_{k} = {\begin{bmatrix} M_{k} & T_{k}^{T} \\ T_{k} & Q_{k} \end{bmatrix} = {\begin{bmatrix} Q_{k + 1} & {\underset{\_}{T}}_{k}^{T} \\ {\underset{\_}{T}}_{k} & {\underset{\_}{M}}_{k} \end{bmatrix}.}}} & (40) \end{matrix}$

[0161] Using this notation, equation (36) of the time steps k and k+1 is included in the following equation. $\begin{matrix} {{{\overset{˘}{Q}}_{k}\begin{bmatrix} 0 \\ K_{k} \end{bmatrix}} = {\begin{bmatrix} \alpha_{k}^{T} \\ C_{k}^{T} \end{bmatrix} = {{\overset{˘}{C}}_{h}^{T} + \begin{bmatrix} {\alpha_{k}^{T} - c_{k}^{T}} \\ 0 \end{bmatrix}}}} & (41) \\ {{{\overset{˘}{Q}}_{k}\begin{bmatrix} K_{k + 1} \\ 0 \end{bmatrix}} = {\begin{bmatrix} C_{k - 1}^{T} \\ \beta_{k}^{T} \end{bmatrix} = {{\overset{˘}{C}}_{k}^{T} + \begin{bmatrix} 0 \\ {\beta_{k}^{T} - c_{k - N}^{T}} \end{bmatrix}}}} & (42) \end{matrix}$

[0162] where,

α_(k) ^(T) =T _(k) ^(T) K _(k) εR ^(1×2), β_(k) ^(T) =T _(k) K _(k+1) ε R ^(1×2).

[0163] Based on the notation, it is more convenient to obtain K^(U) _(k)εR^((N+1)×2), which satisfies the following equation, than to obtain K_(k) directly.

{haeck over (Q)} _(k) {haeck over (K)} _(k) ={haeck over (C)} _(k) ^(T)  (43)

[0164] where,

{haeck over (K)} _(k) =[k _(k+1) ^(T) K _(k) ^(T)]^(T) =[K _(k+1) ^(T) k _(k−N) ^(T)]^(T)  (44)

[0165] To this end, K^(U) _(k)εR^((N+1)×2) can be expressed as shown in equation (46) by using equation (45), obtained from equation (41). $\begin{matrix} {{\overset{˘}{C}}_{k}^{T} = {{{\overset{˘}{Q}}_{k}\begin{bmatrix} 0 \\ K_{k} \end{bmatrix}} - \begin{bmatrix} {\alpha_{k}^{T} - c_{k}^{T}} \\ 0 \end{bmatrix}}} & (45) \\ \begin{matrix} {{\overset{˘}{K}}_{k} = {\begin{bmatrix} m_{k} \\ \mu_{k} \end{bmatrix} = {{{\overset{˘}{Q}}_{k}^{- 1}{\overset{˘}{C}}_{k}^{T}} = {\begin{bmatrix} 0 \\ K_{k} \end{bmatrix} - {{\overset{˘}{Q}}_{k}^{- 1}\begin{bmatrix} {\alpha_{k}^{T} - c_{k}^{T}} \\ 0 \end{bmatrix}}}}}} \\ {= {\begin{bmatrix} 0 \\ K_{k} \end{bmatrix} - {\begin{bmatrix} S_{k}^{- 1} \\ {A_{k}S_{k}^{- 1}} \end{bmatrix}\left\lbrack {\alpha_{k}^{T} - c_{k}^{T}} \right\rbrack}}} \end{matrix} & (46) \end{matrix}$

[0166] Here, K^(U) _(k) is divided into m_(k)εR^(N×2) and μ_(k)εR^(1×2). Also note that α^(T) _(k)−C^(T) _(k)=−(C_(k) ^(T)+A_(k) ^(T)C_(k) ^(T)). Further, assuming that Q^(U) _(k) has an inverse matrix, auxiliary variables A_(k)εR^(N×1) and S_(k)εR satisfy the following equation. $\begin{matrix} {{{\overset{˘}{Q}}_{k}\begin{bmatrix} 1 \\ A_{k} \end{bmatrix}} = {\begin{bmatrix} S_{k} \\ 0 \end{bmatrix}\left( {{\begin{bmatrix} 1 \\ A_{k} \end{bmatrix}S_{k}^{- 1}} = {{\overset{˘}{Q}}_{k}^{- 1}\begin{bmatrix} 1 \\ 0 \end{bmatrix}}} \right)}} & (47) \end{matrix}$

[0167] where, the bottom block of the above equation means T_(k)+Q_(k)A_(k)=0 or T_(k)T=−A_(k) ^(T)Q_(k) ^(T).

[0168] Next, auxiliary variables B_(k)εR^(N×1) and F_(k)εR such as those shown in the following equation (48) are introduced to delete μ_(k) in equation (46) without affecting the top block of C^(T) _(k). Further, subtracting B^(U) _(k)F_(k) ⁻¹μ_(k) from K^(U) _(k) in equation (46) provides equation (49). $\begin{matrix} {{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k}} = {{{\overset{˘}{Q}}_{k}\begin{bmatrix} B_{k} \\ F_{k} \end{bmatrix}} = {\begin{bmatrix} 0 \\ 1 \end{bmatrix}\left( {{\overset{˘}{B}}_{k} = \begin{bmatrix} B_{k} \\ F_{k} \end{bmatrix}} \right)}}} & (48) \\ \begin{matrix} {{{\overset{˘}{K}}_{k} - {{\overset{˘}{B}}_{k}F_{k}^{- 1}\mu_{k}}} = {\begin{bmatrix} m_{k} \\ \mu_{k} \end{bmatrix} - {\begin{bmatrix} {B_{k}F_{k}^{- 1}} \\ 1 \end{bmatrix}\mu_{k}}}} \\ {= \begin{bmatrix} {m_{k} - {B_{k}F_{k}^{- 1}\mu_{k}}} \\ 0 \end{bmatrix}} \end{matrix} & (49) \end{matrix}$

[0169] Then, the left-hand side of equation (49) is multiplied by Q^(U) _(k) from the left to obtain the following equation. $\begin{matrix} \begin{matrix} {{{\overset{˘}{Q}}_{k}\left( {{\overset{˘}{K}}_{k} - {{\overset{˘}{B}}_{k}F_{k}^{- 1}\mu_{k}}} \right)} = {{{{\overset{˘}{Q}}_{k}{\overset{˘}{K}}_{k}} - {{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k}F_{k}^{- 1}\mu_{k}}} = {{\overset{˘}{C}}_{k}^{T} - {\begin{bmatrix} 0 \\ 1 \end{bmatrix}F_{k}^{- 1}\mu_{k}}}}} \\ {= {{\overset{˘}{C}}_{k}^{T} - \begin{bmatrix} 0 \\ {F_{k}^{- 1}\mu_{k}} \end{bmatrix}}} \end{matrix} & (50) \end{matrix}$

[0170] Equation (49) is substituted for the left-hand side of the above equation. Then, equation (43) is expressed as follows: $\begin{matrix} {{{{\overset{˘}{Q}}_{k}\left( {{\overset{˘}{K}}_{k} - {{\overset{˘}{B}}_{k}F_{k}^{- 1}\mu_{k}}} \right)} = {{\overset{˘}{C}}_{k}^{T} - \begin{bmatrix} 0 \\ {F_{k}^{- 1}\mu_{k}} \end{bmatrix}}},{{\begin{bmatrix} Q_{k + 1} & {\underset{\_}{T}}_{k}^{T} \\ {\underset{\_}{T}}_{k} & {\underset{\_}{M}}_{k} \end{bmatrix}\begin{bmatrix} {m_{k} - {B_{k}F_{k}^{- 1}\mu_{k}}} \\ 0 \end{bmatrix}} = {\begin{bmatrix} C_{k + 1}^{T} \\ c_{k - N}^{T} \end{bmatrix} + \begin{bmatrix} 0 \\ {{- F_{k}^{- 1}}\mu_{k}} \end{bmatrix}}}} & (51) \end{matrix}$

[0171] This is the same form as equation (42). The following equation (52) can be obtained from the top block of equation (51).

Q _(k+1)(m _(k) −B _(k) F _(k) ⁻¹ μ _(k))=C _(k+1) ^(T)  (52)

[0172] Equations (36) and (52) are compared to obtain the update equation of the gain matrix K_(k).

[0173] (Lemma 4)

[0174] The auxiliary variables A_(k) and S_(k) can be obtained as follows:

A _(k) =A _(k−1) −K _(k) W _(k)[c_(k) +C _(k) A _(k−1) ]εR ^(N×1)  (53)

S _(k) =ρS _(k−1) +[c _(k) ^(T) +A _(k) ^(T) C _(k) ^(T) ]W _(k) [c _(k) +C _(k) A _(k−1) ]ε R  (54)

[0175] where, A⁻¹=0, and S⁻¹=1/ε₀.

[0176] (Demonstration) By using the following equation (55) of A_(k) and S_(k) and equation (39), equation (56) is obtained. $\begin{matrix} {{{\overset{˘}{Q}}_{k - 1}\begin{bmatrix} 1 \\ A_{k - 1} \end{bmatrix}} = \begin{bmatrix} S_{k - 1} \\ 0 \end{bmatrix}} & (55) \\ \begin{matrix} {{{\overset{˘}{Q}}_{k}\begin{bmatrix} 1 \\ A_{k - 1} \end{bmatrix}} = {{\rho \quad {{\overset{˘}{Q}}_{k - 1}\begin{bmatrix} 1 \\ A_{k - 1} \end{bmatrix}}} + {{\overset{˘}{C}}_{k}^{T}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}} \\ {= {\begin{bmatrix} {\rho \quad S_{k - 1}} \\ 0 \end{bmatrix} + {\begin{bmatrix} c_{k}^{T} \\ C_{k}^{T} \end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}}} \end{matrix} & (56) \end{matrix}$

[0177] On the other hand, the following equation is obtained by multiplying both sides of equation (41) by W_(k) [C _(k)+C_(k)A_(k−1)]. $\begin{matrix} {{{{\overset{˘}{Q}}_{k}\begin{bmatrix} 0 \\ K_{k} \end{bmatrix}}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}} = {\begin{bmatrix} \alpha_{k} \\ C_{k}^{T} \end{bmatrix}{{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}.}}} & (57) \end{matrix}$

[0178] By subtracting equation (57) from equation (56), the following equation (58) is formed. $\begin{matrix} {{{{{{{{{{{{{{\overset{˘}{Q}}_{k}\left\lbrack {\begin{bmatrix} 1 \\ A_{k - 1} \end{bmatrix} - {\begin{bmatrix} 0 \\ K_{k} \end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}} \right\rbrack} =}\quad}\begin{bmatrix} {\rho \quad S_{k - 1}} \\ 0 \end{bmatrix}} +}\quad}\begin{bmatrix} c_{k}^{T} \\ C_{k}^{T} \end{bmatrix}}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}} - {\begin{bmatrix} \alpha_{k}^{T} \\ C_{k}^{T} \end{bmatrix}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}},{{{\overset{˘}{Q}}_{k}\begin{bmatrix} 1 \\ {A_{k - 1} - {K_{k}{W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}} \end{bmatrix}} =}}\quad}\quad}{\quad\begin{bmatrix} {{\rho \quad S_{k - 1}} + {\left\lbrack {c_{k}^{T} - \alpha_{k}^{T}} \right\rbrack {W_{k}\left\lbrack {c_{k} + {C_{k}A_{k - 1}}} \right\rbrack}}} \\ 0 \end{bmatrix}}} & (58) \end{matrix}$

[0179] This equation is compared with equation (47). Since α_(k) ^(T)=T_(k) ^(T)K_(k)=−A_(k) ^(T)C_(k) ^(T), equations (53) and (54) are obtained.

[0180] (Lemma 5)

[0181] The auxiliary variable D_(k)=B_(k)F_(k) ⁻¹ is obtained by the following equation (59). F_(k) is updated by the following equation (60).

D _(k) =[D _(k−1) −m _(k) W _(k)η_(k][1−μ) _(k)W_(k)η_(k)]⁻¹ ε R^(N×1)  (59)

F _(k) =F _(k−1)[1−μ_(k) W _(k)η_(k)]/ρεR  (60)

[0182] where, η_(k)=C^(U) _(k)D^(U) _(k−1)=c_(k−N)+C_(k+1)D_(k−1), D⁻¹=0, and F⁻¹=0.

[0183] (Demonstration) In order to update B_(k) and F_(k), equation (62) is formed by using equation (61). $\begin{matrix} {{{\overset{˘}{Q}}_{k - 1}{\overset{˘}{B}}_{k - 1}} = {{{\overset{˘}{Q}}_{k - 1}\begin{bmatrix} B_{k - 1} \\ F_{k - 1} \end{bmatrix}} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}}} & (61) \\ \begin{matrix} {{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} = {{\rho \quad {\overset{˘}{Q}}_{k - 1}{\overset{˘}{B}}_{k - 1}} + {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}}} \\ {= {{\rho \begin{bmatrix} 0 \\ 1 \end{bmatrix}} + {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}}} \end{matrix} & (62) \end{matrix}$

[0184] In order to modify the above equation so as to have the same form as equation (61), C^(U) _(k) ^(T)W_(k)C^(U) _(k)B_(k−1) ^(U) is subtracted from equation (62) to obtain the following equation. $\begin{matrix} {{{{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} - {{\overset{˘}{C}}_{k}^{T}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}} = {{{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k - 1}} - {{\overset{˘}{Q}}_{k}{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}}} = {\rho \begin{bmatrix} 0 \\ 1 \end{bmatrix}}}},{{{\overset{˘}{Q}}_{k}\left\lbrack {{\overset{˘}{B}}_{k - 1}\overset{˘}{K}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}} \right\rbrack} = {\rho \begin{bmatrix} 0 \\ 1 \end{bmatrix}}}} & (63) \end{matrix}$

[0185] Comparing the above last equation with equation (48) yields a recursive equation for B^(U) _(k).

{haeck over (B)} _(k)=({haeck over (B)} _(k−1) −K _(k) W _(k) {haeck over (B)} _(k−1))/ρ  (64)

[0186] $\begin{matrix} {{D_{k} = {B_{k}F_{k}^{- 1}}},{{\overset{˘}{D}}_{k} = {{{\overset{˘}{B}}_{k}F_{k}^{- 1}} = \begin{bmatrix} D_{k} \\ 1 \end{bmatrix}}}} & (65) \end{matrix}$

[0187] B_(k) and F_(k) are updated by this equation.

[0188] Since they appear only for B^(U) _(k) and D_(k)=B_(k)F_(k) ⁻¹εR^(N×1), however, it is more convenient to express equations (48) and (64) to the following equation (65). The matrix D_(k) satisfies the following equation (66). $\begin{matrix} {{{{\overset{˘}{Q}}_{k}{\overset{˘}{D}}_{k}} = {{{\overset{˘}{Q}}_{k}{\overset{˘}{B}}_{k}F_{k}^{- 1}} = {\begin{bmatrix} 0 \\ 1 \end{bmatrix}F_{k}^{- 1}}}},\quad {{{\overset{˘}{Q}}_{k}\begin{bmatrix} D_{k} \\ 1 \end{bmatrix}} = \begin{bmatrix} 0 \\ F_{k}^{- 1} \end{bmatrix}}} & (66) \\ \begin{matrix} {{{\overset{˘}{Q}}_{k}\left\lbrack {{{\overset{˘}{B}}_{k - 1}F_{k - 1}^{- 1}} - {{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{B}}_{k - 1}F_{k - 1}^{- 1}}} \right\rbrack} = {{\overset{˘}{Q}}_{k}\left\lbrack {{\overset{˘}{D}}_{k - 1} - {{\overset{˘}{K}}_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack}} \\ {= \begin{bmatrix} 0 \\ {\rho \quad F_{k - 1}^{- 1}} \end{bmatrix}} \end{matrix} & (67) \end{matrix}$

[0189] Next, equation (63) is multiplied by F_(k−1) ⁻¹ to obtain the following equation (67), and is further expressed by the following equation (68) when D^(U) _(k−1)=B^(U) _(k−1)F_(k−1) ⁻¹ is used. $\begin{matrix} {{{{\overset{˘}{Q}}_{k}\left\lbrack {{\overset{˘}{D}}_{k - 1} - {\begin{bmatrix} m_{k} \\ \mu_{k} \end{bmatrix}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \right\rbrack} = \begin{bmatrix} 0 \\ {\rho \quad F_{k - 1}^{- 1}} \end{bmatrix}},{{{\overset{˘}{Q}}_{k}\begin{bmatrix} {D_{k - 1} - {m_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \\ {1 - {\mu_{k}W_{k}{\overset{˘}{C}}_{k}{\overset{˘}{D}}_{k - 1}}} \end{bmatrix}} = \begin{bmatrix} 0 \\ {\rho \quad F_{k - 1}^{- 1}} \end{bmatrix}}} & (68) \end{matrix}$

[0190] Therefore, the following equation is obtained when equation (68) is multiplied by [1−μ_(k)W_(k)C^(U) _(k)D^(U) _(k−1)]⁻¹. $\begin{matrix} \begin{matrix} {{\overset{\bigvee}{Q}}_{k}\begin{bmatrix} {\left\lbrack {D_{k - 1} - {m_{k}W_{k}{\overset{\bigvee}{C}}_{k}{\overset{\bigvee}{D}}_{k - 1}}} \right\rbrack \left\lbrack {1 - {\mu_{k}W_{k}{\overset{\bigvee}{C}}_{k}{\overset{\bigvee}{D}}_{k - 1}}} \right\rbrack}^{- 1} \\ 1 \end{bmatrix}} \\ {\quad {= \begin{bmatrix} 0 \\ {\rho \quad {F_{k - 1}^{- 1}\left\lbrack {1 - {\mu_{k}W_{k}{\overset{\bigvee}{C}}_{k}{\overset{\bigvee}{D}}_{k - 1}}} \right\rbrack}^{- 1}} \end{bmatrix}}} \end{matrix} & \quad \end{matrix}$

[0191] By comparing this equation with equation (66), an update equation for D_(k) and F_(k) is finally obtained.

[0192] (Demonstration of Lemma 6 (Existence Condition Suitable for Fast Processing))

[0193] As described above, the existence of the fast H_(∞) filter can be checked with the computational complexity of O(N) by using the existence condition of equations (22) and (23). A demonstration thereof will be shown below.

[0194] When the characteristic equation of a 2×2 matrix R_(e, k) shown in the following equation (69) is solved, a eigenvalue λ_(i) of R_(e, k) is obtained by the following equation (70). $\begin{matrix} \begin{matrix} {{{{\lambda \quad I} - R_{e,k}}} = {\begin{matrix} {\lambda - \left( {\rho + {H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}}} \right)} & {{- H_{k}}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}} \\ {{- H_{k}}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}} & {\lambda - \left( {{{- \rho}\quad \gamma_{f}^{2}} + {H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}}} \right)} \end{matrix}}} \\ {= {{\lambda^{2} - {\left( {{2H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}} + {\rho \quad \varrho}} \right)\lambda} - {\rho^{2}\gamma_{f}^{2}} + {\rho \quad \varrho \quad H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}}} = 0}} \end{matrix} & (69) \\ {\lambda_{i} = \frac{\Phi \pm \sqrt{\Phi^{2} - {4\rho \quad \varrho \quad H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}} + {4\rho^{2}\gamma_{f}^{2}}}}{2}} & (70) \end{matrix}$

[0195] where, Φ=2H_(k){circumflex over (Σ)}_(k|k−1)H_(k) ^(T)+ρe, e=1−γ_(f) ²

[0196] If the following expression (71) is satisfied, one of the two eigenvalues of the matrix R_(e, k) is positive and the other is negative, and the matrixes R_(k) and R_(e, k) have the same inertia. Therefore, the existence condition of equation (22) is obtained by using the following equation (72). Here, the calculation of H_(k)K{tilde over ()}_(k) requires the same number of multiplications as O(N).

−4ρeH _(k){circumflex over (Σ)}_(k|k−1) H _(k) ^(T)+4ρ²γ_(f) ²>0  (71 )

[0197] $\begin{matrix} {{H_{k}{{\hat{\sum}}_{k|{k - 1}}\quad H_{k}^{T}}} = \frac{H_{k}{\overset{\sim}{K}}_{b}}{1 - {H_{k}K_{b}}}} & (72) \end{matrix}$

[0198] Industrial Applicability

[0199] According to the present invention, as described above, the fast real-time identification of time-invariant and time-variant systems can be implemented by using the fast algorithm (fast H_(∞) filtering algorithm) for the modified H_(∞) filters developed based on the new H_(∞) evaluation criterion. In addition, according to the present invention, the present algorithm includes, as a particular case, the fast Kalman filtering algorithm, and a term corresponding to the covariance of system noise which is dominant in the tracking performance of a time-varying system can be theoretically determined. Further, according to the present invention, a fast time-varying system identification method can be provided, which is very effective particularly when a system (impulse response) is varied discontinuously with time, such as an echo canceller for a time-varying system which varies extremely as sudden line switching. Furthermore, according to the present invention, a system identification method can be provided, which is applicable to echo cancellers in communication systems and acoustic systems, sound-field reproduction, and noise control. 

1. A system identification method for performing a fast real-time identification of a time-invariant or time-variant system, wherein an H_(∞) filter equation expressed by the following equation is used, {circumflex over (x)} _(k+1)|k+1 =x _(k|k) +K _(s, k+1)(γ_(k+1) H _(k+1) x _(k|k)) where, {circumflex over (x)}_(k|k): The estimate of state x_(k) at time k, obtained by using observation signals y_(o) to y_(k) y_(k): Observation signal K_(s, k+1): Filter gain H_(k): Observation matrix and a filtering algorithm robust against disturbance is formed by setting, as an H_(∞) evaluation criterion, a maximum energy gain from disturbance weighted by a weight (ρ) of the evaluation function to a filter error to be smaller than a predetermined upper limit (γ_(f) ²)
 2. A system identification method according to claim 1, wherein, as the H_(∞) evaluation criterion, the maximum value of {(value indicating the filter error/weight (ρ) of the evaluation function)/[(value indicating an initial state)+(value indicating system noise)+(value indicating observation noise/weight (ρ) of the evaluation function)]} is set to be smaller than the predetermined upper limit (γ_(f) ²).
 3. A system identification method according to claim 1 or 2, wherein an output signal is further obtained from the state estimate x{circumflex over ( )}_(k|k) at time k by the following equation. {circumflex over (z)} _(k|k) =H _(k){circumflex over (x)}_(k|k) z_(k): Output signal
 4. A system identification method according to any one of claims 1 to 3, wherein the following equation shows the relationship between the weight (ρ) of the evaluation function and the predetermined upper limit (γ_(f) ²) 0<ρ=1−γ_(f) ⁻²≦1γ_(f)>1
 5. A system identification method according to any one of claims 1 to 4, wherein the filter gain K_(s, k) is given by the following relational equation by using a gain matrix K_(k). K_(s,k) =G _(k) ⁻¹ {tilde over (K)} _(k) , G _(k)=ρ+γ_(f) ⁻² H _(k) {tilde over (K)} _(k) εR where, {tilde over (K)} _(k)(i)=ρK _(k)(i, 1), i=1,2, . . . ,N.
 6. A system identification method according to any one of claims 1 to 5, comprising steps of; setting initial conditions of a recursive equation of the gain matrix K_(k), the auxiliary variables, and the state estimate x{circumflex over ( )}_(k|k), recursively determining the auxiliary variables at time k and obtaining a second gain matrix in which a row including the auxiliary variables is added to the gain matrix K_(k), dividing the second gain matrix and obtaining first and second divisional gain matrixes, obtaining a gain matrix K_(k+1) at time k+1 from an equation including the first and second divisional gain matrixes, and obtaining a filter gain K_(s, k+1) at time k+1 from the relational equation of the gain matrix K_(k) and the filter gain K_(s, k), updating the H_(∞) filter equation according to the obtained filter gain K_(s, k+1), and repeating each of the above steps with the time being put forward.
 7. A system identification method according to any one of claims 1 to 6, wherein the existence of the fast H_(∞) filter is checked by using the following equation as an existence condition suitable for fast processing, with the computational complexity of O(N). −e ^({circumflex over (Ξ)}) _(i)+ργ_(f) ²>0, i=0, . . . ,k where, ${\varrho = {1 - \gamma_{f}^{2}}},{{\hat{\Xi}}_{i} = \frac{H_{i}{\overset{\sim}{K}}_{i}}{1 - {H_{i}{\overset{\sim}{K}}_{i}}}}$


8. A system identification method according to any one of claims 1 to 7, wherein an echo canceller is implemented by applying the H_(∞) filter equation to obtain the state estimate x{circumflex over ( )}_(k|k), producing a quasi echo as in the following equation, and canceling an actual echo by the obtained quasi echo. ${\hat{d}}_{k} = {{H_{k}{\hat{x}}_{k|k}} = {\sum\limits_{i = 0}^{N - 1}\quad {{{\hat{h}}_{i}\lbrack k\rbrack}u_{k - i}}}}$

H _(k)=[u_(k) , . . . ,u _(k−N+1)] where, {circumflex over (d)}_(k) Quasi-echo u_(k) Received signal N Tap number ĥ_(i)[k] Estimate of impulse response of echo path 